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Abstract. In this paper, we propose a distributed algorithm for solving large-scale 



separable convex problems using Lagrangian dual decomposition and the interior-point 



framework. By adding self-concordant barrier terms to the ordinary Lagrangian, we 



prove under mild assumptions that the corresponding family of augmented dual func- 



tions is self-concordant. This makes it possible to efficiently use the Newton method 



for tracing the central path. We show that the new algorithm is globally convergent 



and highly parallelizable and thus it is suitable for solving large-scale separable convex 



problems. 



Keywords. Separable convex optimization, self-concordant functions, interior-point 



methods, augmented Lagrangian decomposition, parallel computations. 
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1 Introduction 



Can self- concordance and interior-point methods be incorporated into a Lagrangian 



dual decomposition framework? This paper presents a decomposition algorithm that 



incorporates the interior-point method into augmented Lagrangian decomposition tech- 



nique for solving large-scale separable convex problems. Separable convex problems, 



i.e. optimization problems with a separable convex objective function but with coupling 



constraints, arise in many fields: networks (communication networks, multicommodity 



network flows) [HIS], process system engineering (e.g. distributed model predictive con- 



trol) [SIS], stochastic programming [S], etc. There has been considerable interest in par- 



allel and distributed computation methods for solving this type of structured optimiza- 



tion problems and many methods have been proposed: dual subgradient methods [SIE], 



alternating direction methods [6l[8], proximal method of multipliers [9], proximal center 



method [1], interior-point based methods [2|l5| [T0] - [n] . etc. 



The methods mentioned above belong to the class of augmented Lagrangian or 
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multiplier methods [6], i.e. they can be viewed as techniques for maximizing an aug- 



mented dual function. For example in the alternating direction method a quadratic 



penalty term is added to the standard Lagrangian to obtain a smooth dual function 



and then using a steepest ascent update for the multipliers. However, the quadratic 



term destroys the separability of the given problem. Moreover, the performance of 



these methods is very sensitive to the variations of their parameters and some rules for 



choosing these parameters were given e.g. in [8l[T5]. In the proximal center method [1] 



we use smoothing techniques in order to obtain a well-behaved Lagrangian, i.e. we add 



a separable strongly convex term to the ordinary Lagrangian. This technique leads to a 



smooth dual function, i.e. with Lipschitz continuous gradient, which preserves separa- 



bility of the problem, the corresponding parameter is selected optimally and moreover 



the multipliers are updated using an optimal gradient based scheme. In p,[5l, [T0] - [ll] 



interior-point methods are proposed for solving special classes of separable convex prob- 



lems with a particular structure of the coupling/local constraints. In those papers the 



Newton direction is used to update the primal variables and/or multipliers obtaining 
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polynomial-time complexity for the proposed algorithms. 



In the present paper we use a similar smoothing technique as in |1] in order to ob- 



tain a well-behaved augmented dual function. Although we relax the coupling con- 



straints using the Lagrangian dual framework as in the main difference here is that 



the smoothing term is a self- concordant barrier, while in [1] the main property of the 



smoothing term was strong convexity. Therefore, using the properties of self-concordant 



functions we show that the augmented dual function becomes under mild assumptions 



also self-concordant. Hence the Newton direction can be used instead of gradient based 



directions as it is done in most of the augmented Lagrangian methods. Furthermore, we 



develop a specialized interior-point method to maximize the augmented dual function 



which takes into account the special structure of our problem. We present a parallel 



algorithm for computing the Newton direction of the dual function and we also prove 



global convergence of the proposed method. 



The main contributions of the paper are the following: 



(i) We consider a more general model for separable convex problems that includes local 
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equality and inequality constraints, and linear coupling constraints which generalizes 
the models in [2l[5l[T0l[Tl] . 

(ii) We derive sufficient conditions for self-concordance of augmented Lagrangian and 
we prove self-concordance for the corresponding family of augmented dual functions. 

(iii) We provide an interior-point based algorithm for solving the dual problem with 
proofs of global convergence and polynomial-time complexity. 

(iv) We propose a practical implementation of the algorithm based on solving approx- 
imately the subproblems and on parallel computations of the Newton directions. 

Note that item (ii) generalizes the results of [5l[T0]. However, the consideration of 
general convex problems with local equality constraints requires new proofs with more 
involved arguments in order to prove self- concordance for the family of dual functions. 

This paper is organized as follows. In Section [2] we formulate the separable convex 
problem followed by a brief description of some of the existing decomposition methods 
for this problem. The main results are given in Sections [3] and HI In Section [3] we show 
that the augmented Lagrangian obtained by adding self-concordant barrier terms to the 
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ordinary Lagrangian forms a self-concordant family of dual functions. Then an interior- 



point Lagrangian decomposition algorithm with polynomial complexity is proposed in 



Section HI The new algorithm makes use of the special structure of our problem so that 



it is highly parallelizable and it can be effectively implemented on parallel processors. 



We conclude the paper with some possible applications. 



Throughout the paper we use the following notations. For a function tp with two 



arguments, scalar parameter t and decision variable x, i.e. ip{t,x), we use " /" to 



denote the partial derivative of ip{t, x) with respect to t and "V" with respect to x: e.g. 



Vifj^tjx) = -Mw^ipit.x). For a function (j), three times differentiable, i.e. (j) G C^(dom{/)), 



V'^0(x)[/ii, /i2, /is] denotes the third differential of at x along directions hi, /12 and /13. 



We use the notation A ^ B ii B — Ais positive semidefinite. We use Da to denote the 



block diagonal matrix having on the main diagonal the matrices Ai, ■ ■ ■ , Ajy. We use 



int(X) to denote the interior of a set X. 
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2 Problem Formulation 



We consider the following general separable convex optimization problem: 



N 

r= ,^niin (1) 

i=l 

N 

s.t. BiXi = &, AiXi = ttiWi = 1 ■ ■ ■ N, (2) 



where fi : M"* — )■ M are convex functions, Xi are closed convex sets, Ai G W' 



Bi G M'"^"-', ai G M'"^ and 6 G M™. For simplicity of the exposition we define the vector 
X := [x{---x]^Y, the function f{x) : = ^^e set X := Hili^i) the 

matrix i? := ■ ■ ■ En] and n := X^^;^ ^ij. 



Remark 2.1. (i) Note that we do not assume strict /strong convexity of any function fi. 



(ii) Coupling inequality constraints ^^^^ BiXi < b can be included in this framework 
by adding a slack variable xn+i- J2iLi ^i^i + ^^a^+i = b, i.e. Bn+i = I and X^+i = W]^. 



Let (■,■) denote the Euclidian inner product on M". By forming the Lagrangian 
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corresponding only to the coupling linear constraints (with the multipliers A G M™), i.e. 



Lo{x,\) = f{x) + {\,Bx-b), 



we can define the standard dual function 



(io(A) = min{Lo(x, A) : G Xj, AiXi = Vi = 1 ■ ■ ■ A^}. 

X 



Since Lq preserves the separability of our problem we can use the dual decomposition 



method [6] by solving in parallel N minimization problems and then updating the 



multiphers in some fashion. Note that the dual function do is concave but, in general do 



is not differentiable (e.g. when / is not strictly convex). Therefore, for maximizing do 



one has to use involved nonsmooth optimization techniques [6l[7] . From duality theory 



one knows that if {x*, A*) is a saddle point for the Lagrangian Lq, then under appropriate 



conditions (constraint qualification), x* is an optimal solution for the primal ([I])-([2]) 
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and A* is an associated dual optimal multiplier for the dual problem: maxAGR™ doW- 



In order to obtain a smooth dual function we need to use smoothing techniques 



applied to the ordinary Lagrangian Lq. One approach is the augmented Lagrangian 



obtained e.g. by adding a quadratic penalty term to the Lagrangian Lq: — a|p. In 



the alternating direction method [6],[8] the minimization of the augmented Lagrangian is 



performed by alternating minimization in a Gauss-Seidel fashion followed by a steepest 



ascent update for the multipliers. 



In [1] we proposed the proximal center method in which we added to the standard 



Lagrangian a smoothing term t Yli=i 9Xi{xi), where each function is strongly convex 



and depends on the set Xj so that the augmented Lagrangian takes the following form: 



N 

Lr^(x, A) = Y^iMx,) + tgxAx^)] + (A, Bx-b). 

i=l 



Therefore, the augmented Lagrangian L^^°^ is strongly convex, preserves separability 
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of the problem like Lq and the associated augmented dual function 



df^°^{X) = min{L^™'^(a;, X) : Xi E Xi, AiXi = V« = 1 ■ ■ ■ A^} 



is differentiable and has also a Lipschitz continuous gradient. In |1] an accelerated 



gradient based method is used to maximize the augmented dual function while the 



corresponding minimization problems are solved in parallel. Moreover, the smoothing 



parameter t is selected optimally. 



Note that the methods discussed above use only the gradient directions of the aug- 



mented dual function in order to update the multipliers. Therefore, in the absence 



of more conservative assumptions like strong convexity, the global convergence rate of 



these methods is slow, in general sub-linear. In this paper we propose to smoothen 



the Lagrangian by adding instead of strongly convex terms gx,, self-concordant barrier 
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terms (j)Xi associated with the sets Xi, in order to obtain the self- concordant Lagrangian: 



AT 

Lr{x,X) = + {X,Bx - h). (3) 

i=l 



In the next section we show, using the theory of self-concordant barrier functions [Mp!7] . 



that for a relatively large class of convex functions (see also Section [5]), we can obtain 



a self-concordant augmented dual function: 



(F\t, A) = min{Lf (x, A) : G int(Xi), A^Xi = Vi = 1 ■ ■ ■ iV}. (4) 



This opens the possibility of deriving an interior-point method using Newton di- 



rections for updating the multipliers to speed up the convergence rate of the proposed 



algorithm. 
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3 Sufficient Conditions for Self-Concordance of the 
Augmented Dual Function 

In this section we derive sufficient conditions under which the family of augmented 
dual functions is self-concordant. A key property that allows to prove polynomial 
convergence for barrier type methods is the property of self-concordance (see Definition 
2.1.1 in [16]): 

Definition 3.1. A closed convex function with open convex domain C M" is 
called M(i)- self-concordant, where > 0, if (/) is three times continuously different iable 
on and if for all x G X^ and h G we have 



VV(x)[/i, h, h] < M^{h^W^(t){x)h) ' . (5) 



A function is called N^- self -concordant harrier for its domain X^ if is 2-self- 
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concordant function and for all x G and h G M" we have 



(V0(a;), hf < h^W^(P{x)h. (6) 



Note that (E]) is equivalent to (see [TB], pp. 14): 



3 
i=l 



Moreover, if Hessian V^0(a;) is positive definite, then the inequality ([6]) is equivalent to 



V0(x)^VV(a;)"V0(x) < N^. (8) 



Next lemma provides some basic properties of self-concordant functions: 



Proposition 3.1. ( [TSf . pp. 15) Let (j) be an M^-self-concordant function such that 



its domain does not contain straight lines (i.e. sets of the form {x + au : a E M.}, 



where x G X^ and u ^ 0). Then, the Hessian V^0(x) is positive definite for all x G 
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and (j) is a barrier function for X^. 



□ 



Note that a self-concordant function which is also a barrier for its domain is called 



strongly self- concordant. The next lemma gives some helpful composition rules for self- 



concordant functions. 



Lemma 3.1. (i) [TSf Any linear or convex quadratic function is 0- self- concordant. 



(a) fTSf Let (pi be Mi-self concordant and let pi > 0, i = 1,2. Then the function 



Pi(pi + P202 is also M-self concordant, where M = ma.x{Mi/ ^/p^, M2/ a/P2}- 

(Hi) Let Xf)ox = nr=i[^«' such that li < Ui and tp G C^{int{Xiyox)) be convex. If there 



exists /3 > such that for all x G int{Xhox) o.nd h eW^ the following inequality holds 



then il)t{x) = ip{x) — i ^"=ilog(ui — Xi){xi — U) is 2(1 + (3)/\/t-self concordant. 



Proof, (i) and (ii) can be found in [16], pp. 13. 



n 




(9) 
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(iii) Denote (j)hox{x) = — X]r=i ^'^si'^i ~ ~ h)- Note that 



i=l 



V'0box(x)[/i, h,h]=2Y^ hy{u, - Xif - h^Jixi - kf 



i=l 



and using Cauchy-Schwarz inequality it follows that 0box is 2-self-concordant func- 



tion on int(Xbox)- Let us denote 



c = \/h'^'V'^ip{x)h and d 



\ 



Using hypothesis (|9]) and 2-self-concordance of </)box we have the following inequahties: 



\V'^^t{x)\h,hM < \V^ij{x)[h,h,h]\+t\V^^i,oAx)[h,h,h]\ <(3c'^d + 2td^ 
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With some computations we can observe that 



and since 



h^W'^^Pt{x)h = + td^, 



the proof is complete. □ 
Note that condition ([9]) is similar to if) is (i-compatihle with 0box on Xbox, defined in |16j . 
The following assumptions will be valid throughout this section: 

Assumption 3.1. We consider a given compact convex set X with nonempty interior 
and an associated A^<^-self-concordant barrier for X (whenever X = Xbox we consider 
= 0box)- Given a function / G C^(int(X)), we also assume that it satisfies one 



of the properties (i)-(iii) of Lemma I3.H i.e. / is either linear or convex quadratic or 
Mj-self-concordant or X is a box and / satisfies condition Let A G < n, 
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and B G R™^" be so that the matrix 



Ax = a} nint(X)} 



A 



B 



has full row rank and the set {{x G 



We analyze the following prototype minimization problem: 



min{/(x) + t(j){x) + (A, Bx) : x G int(X), Ax = a}. 



(10) 



Let us define the dual convex function: 



d(t, A) := maxj— /(x) — — (A, Bx) : x G int(X), Ax = a}. 



Boundedness of the set X and self-concordance property of the function f + t(j) (which 



follow from the assumptions mentioned above) guarantee existence and uniqueness of 



the maximizer x(t, A) of ffTOl) . Therefore, we can consistently define the maximizer 



x{t, A) of ffTUl) and the dual convex function d{t, A) for every t > and A G 



In the following four lemmas we derive the main properties of the family of aug- 



mented dual functions {d{t, ■)}t>o- We start with a linear algebra result: 



Lemma 3.2. Let A e RP^",p < n, and B G R™-^" he two matrices and U be the matrix 



whose columns form a basis of the null space of A. Then the matrix 



row rank if and only if BU and A have full row rank. 



A 



B 



has full 



Proof. Assume that 



A 



B 



has full row rank. Then A has full row rank. It remains to 



show that BU has full row rank. Assume that this is not the case then there exists a 



vector X G M™, x 7^ such that x^BU = 0. Since the columns of U span the null space 



of A which is orthogonal on the image space of A'^ , it follows that x^B belongs to the 



image space of A , i.e. there exists some y eMP such that x B = y A. But from the 



fact that 



A 



B 



has full row rank we must have x = and y = which contradicts 



our assumption on x. 



If BU and A have full row rank, it follows immediately that B must have full row rank. 
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Assume that 



A 



B 



does not have full row rank. Since A has full row rank, then there 



exist some y G M"* and x G x ^ 0, such that 



y^A + x^B = 0. 



It follows also that 



y^AU + x^BU = 0, i.e. x' BU = 



and thus x = which is a contradiction. 



□ 



Lemma 3.3. // Assumption \3.1\ holds, then for any t > the function d{t, ■) is Mt- 



self- concordant, where Mt is either 2/ \/t or max 



Proof. Since / is assumed to be either linear or convex quadratic or Mj-self-concordant 



or X is a box and / satisfies condition Q it follows from Lemma 13.11 that / + t0 is 



also Mt-self concordant (where Mt is either 2/-\/t or max{Mj, 2/A/t} or 2(1 + f3)/^/t. 



20 



respectively) and with positive definite Hessian (according to our assumptions and 



Proposition 13. ip . Moreover, f + t(f) is strongly self-concordant since is a barrier 



function for X. Since A has full row rank and p < n, then there exists some vectors 



Ui,i = 1 ■ ■ -n — p, that form a basis of the null space of this matrix. Let U be the 



matrix having as columns the vectors Ui and xq a particular solution of Ax = a. Then, 



for a fixed t, the feasible set of fITU]) can be described as 



Q = {y eR'' : Xo + Uy e int(X)}, 



which is an open convex set. Using that self-concordance is affine invariant it follows 



that the functions f{y) = /(xq + Uy), (f){y) = 0(xo + Uy) have the same properties as 



the functions /, 0, respectively, that f + t(j) is also M^-self concordant and that 



d{t, A) = max[-/(y) - t^{y) - (A, B{xo + Uy))]. 
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From our assumptions and Proposition 13 . 1 1 it follows that the Hessian of </> and (f) are 



positive definite. Since / is convex it follows that the Hessian oi f + t(f) is also positive 



definite and thus invertible. Let 



F{t, A) = max[(A, y) - f{y) - t^{y)] 



be the Legendre transformation of / + t(f). In view of known properties of the Legendre 



transformation, it follows that if / + t0 is convex on X from such that its Hessian 



is positive definite, then F{t,-) has the same properties on its domain {A G M™ : 
(A, y) — f{y) — t(j){y) bounded above on Q}. Moreover, from Theorem 2.4.1 in it 



follows that F(t, ■) is also Mf-self-concordant on its domain. Note that 



d{t, A) = (A, -Bxo) + -{BUfX). 
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Since 



A 



B 



has full row rank, then from Lemma 13.21 BU has full row rank. Moreover, 



since V^-F(t, ■) is positive definite and 



V'^d{t,\) = BUV^F{t,-{BUfX){BUf 



it follows that V'^d{t, ■) is positive definite on its domain 



Xd(t,-) ■■= {A G R"" : -f{y) - t(p{y) - (A, B{xq + Uy)) bounded above on Q}. 



Moreover, since self- concordance is afiine invariant it follows that d{t, ■) is also Mr self- 



concordant on the domain ^^(4,.). 



□ 



Lemma 3.4. Under As sumption ] 3. 1\ the inequality \ (yd'(t. X).h) \ < {2^t/ ott)\Jh'^V^d{t, X)h 



holds true for each t > and A, /i G M™, where = {Mt/2)^J N^/t and at = Mt. 



Proof. From Lemma [3.31 we know that d{t,-) is with positive definite Hessian. By 
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virtue of the barrier <^ for the set X the optimal solution A) of (ITOl) satisfies x{t, A) G 
int(X) and so the first-order optimality conditions for optimization problem (fTOj) are: 
there exists u{t, A) G M*' such that 



V/(x(t, A)) + tV0(x(t, A)) + B'^X + A) = and Ax{t, A) = a. (11) 



First we determine the formula for the Hessian. It follows immediately from ( ITTj) that 



Vd(t, A) = -Bx{t, A) and V^rf(t, A) = -B\/x{t, A). 



Let us introduce the following notation: 



H{t, A) := V^fix{t, A)) + tV^0(x(t, A)). 



For simplicity, we drop the dependence of all the functions on x(t, A) and (t, A). Differ- 
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entiating ( fTTl) with respect to A we arrive at the following system in Vx and Vz/: 



V2/ + tV2(/. 




Va; 






A 




Vz/ 








Since H is positive definite and according to our assumption A is full row rank, it 



follows that the system matrix is invertible. Using the well-known formula for inversion 



of partitioned matrices we find that: 



W^d = B[H-^ - H~'A^{AH-'A^y'AH^']B'^. (12) 



Differentiating the first part of (fTTl) with respect to t and using the same procedure as 



before we arrive at a similar system as above in the unknowns x' and u' . We find that 



x' = -[H-^ - H-^A^{AH-^A^)-^AH-^]S/^ and Wd' = -Bx . 
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We also introduce the following notation: F := H"^ {AH~^ A^)~^ AH^^ and G : = 
H^^ — F, which are positive semidefinite. Using a similar reasoning as in [5] and 
Cauchy-Schwarz inequality we obtain: 



Since 

G = H~^-F <H~^ = (VV + tVV)"^ ^ 1A(V^0)" 



and using ([8]) it follows that 



\{Vd',h)\ < .JhT{S/H)h^[Njt. 



□ 



Lemma 3.5. Under Assumption \3. 1\ the inequality |(V^c/'(t, \)h, h)\ < 2r]t h^'V'^d{t, \)h 
holds true for each t>0 andX,he M™, where 7]t = {Mt/2)^ N^/t + {l/2t) . 
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Proof. We recall that H{t, A) = V^f{x{t, A)) + tV^(j){x{t, A)). Therefore 

h^H'it, X)h = {Vy{x{t, A)) + tV^(j){x{t, X)))[x'{t, A), h, h] + h'^V^(j){x{t, X))h, 
h^{H-\t, \)yh = -h^H-\t, X)H'{t, X)H~\t, X)h. 

We again drop the dependence on (t, A) and after some straightforward algebra 
putations we arrive at the following expression: 

{V^d'h, h) = -h^B{H-^ - F)H'{H~^ - F)B^h. 

Let us denote with 

u := {H-^ - F)B^h. 
Taking into account the expression of H' derived above we obtain: 

\{V^d'h,h)\ = lu'^H'ul = liV'^f + tV^(j))[x',u,u] + u^V^(j)u\. 



Using the self- concordance property (JTj) for f + t(j) we obtain that: 



I (VV + iV^0) W, u,u\\< Mtu^iV^f + tV^<P)u^J{x'Y{V^f + tV^<P)x' 

= MtU^Hu^{x'YHx'. 



Moreover, since / is convex, V^/ is positive semidefinite and thus: 



M^V^0M < il/t) u^Hu. 



Combining the last two inequalities we obtain: 



(V^ci'/i, h) I < MtU^Hu^/{x'YHx' + {l/t)u^Hu. (13) 



With some algebra we can check that the following identity holds: FH{H ^ — F) = 0. 



Based on this identity we can compute u^Hu and {x')^Hx'. Indeed, 
u^Hu = h^B{H^^ - F)H{H~^ - F)B^h 

= h^B{H-^ - F)B^h - h^BFH{H-^ - F)B^h = h^B{H-^ - F)B^h = h^W^d h. 



Similarly, using ^ we obtain 



x'fHx' = Vf{H-^ - F)V(j) < VfH-^V(t) < (l/t)V0^(V^0)- V(/> < N^/t. 



The inequality from lemma follows then by replacing the last two relations in fll3p . □ 



The main result of this section is summarized in the next theorem. 



Theorem 3.1. Under the Assumption \3.1\ {c/(t, A)}j>o is a strongly self-concordant 



family in the sense of Definitioi^ 3.1.1 in IW^ with parameters = Mt,^t = {Mt/2)^JTj^Jt 



an 



dr]t = {Mt/2)^/N^ + (l/2t), where Mt is defined in Lemma{3lE 



^Note that according to Definition 3.1.1 in [16] 7* = 1 and /it = 1 for our case. 
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Proof. Basically, from Definition 3.1.1 in [16] we must check three properties: self- 



concordance of d{t, A) (Lemma 13. 3p and that the first and second order derivative of 



d{t, ■) vary with t at a rate proportional to the derivative itself (Lemmas 13.41 and 13. 5p . 



In conclusion, the Lemmas I3.3H3.5I prove our theorem. □ 



It is known [16] that self-concordant families of functions can be minimized by path- 



following methods in polynomial time. Therefore, this type of family of augmented dual 



functions {d{t, ■)}t>o plays an important role in the algorithm of the next section. 



4 Parallel Implementation of an Interior- Point Based 



Decomposition Method 



In this section we develop an interior-point Lagrangian decomposition method for the 



separable convex problem given by ([I])-(l2]). Our previous Theorem 13.11 is the major 



contribution of our paper since it allows us to effectively utilize the Newton method 



for tracing the trajectory of optimizers of the self-concordant family of augmented dual 
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functions 



4.1 Interior-Point Lagrangian Algorithm 



The following assumptions for optimization problem ([I])-(l2]) will be valid in this section: 



Assumption 4.1. (i) The sets Xi are compact convex sets with nonempty interior and 



')x, are A^j-self-concordant barriers for Xj. 



(ii) Each function fi is either linear or convex quadratic or M/--self-concordant or Xi is 



a box and fi satisfies condition 



(iii) The block-angular matrix 



Oi, Bx = 6} nint(X)} ^ 0. 



B 



has full row rank and the set | {x G : AiXi 



Note that boundedness of the set Xi can be relaxed to Xi does not contain straight 



lines and the set of optimal solutions to problem is bounded. Note also that the 



rank assumption (iii) is not restrictive since we can eliminate the redundant equalities 



(see also Lemma [3.21 for other less restrictive conditions). The constraint qualification 
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condition from Assumption 14.11 (iii) guarantees that strong duality holds for problem 
([I])-© and thus there exists a primal-dual optimal solution {x*, A*). 
Let us introduce the dual function: 

d{t, A) = max{-Lf (x, A) : G int(Xj), AiXi = ai\fi = 1 ■ ■ ■ N} 

X 

N 

= (A, 6) + y^max{-/j(xi) - tcpxA^i) - {\, BiXi) : Xi E int(Xj), AiXi = a,} 

i=l 
N 

= {X,b) + Y,di{t,X). 

1=1 

Note that the function d{t, ■) can be computed in parallel by decomposing the original 
large optimization problem ([T])-([2]) into independent small convex subproblems. 

Lemma 4.1. (i) The family {di{t, ■)}t>o is strongly self-concordant with the parameters 
ai{t) = Mi{t),Ut) = {Mi{t)/2)^Njt and r]i{t) = {M,{t)/2)y^Njt + {l/2t), where 
Mi{t) is either 2/ y/t or max{Mf^,2/ y/t} or 2(1 + /3) / y/t for all i = 1---N. 
(a) The family {d{t, ■)}t>o is strongly self- concordant with parameters a{t) = a/y/t, 
^{t) = ^/t and r]{t) = rj/t, for some fixed positive constants a,^ and t] depending on 
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Proof, (i) is a straightforward consequence of Assumption 14.11 and Theorem 13.11 

(ii) Note that d{t,X) = (A, a) + J2iLi^i{^^ '^)- From Proposition 3.1.1 in [16] we 
have that the sum of strongly self-concordant family of functions is also strongly self- 
concordant family with the parameters: a{t) > max,j{a;j(t)} is a positive continuously 
differentiable function on IR+, C,{t) = a{t) meiXi{2^i(t) / ai{t)} and ri{t) = maxi{r]i{t)} . 
Since ai{t) is either 2/y/i or max{Mf^,2/Vt} or 2{l + (3)/y/i for alH = 1 ■ ■ ■ N, it follows 
that we can always choose a(t) = a/\/t, where a = 2 or a = 2(1 + (3). Similarly, we 
can show that there exists positive constants ^ and r] depending on Ni, Mf- and /3 such 
that ^{t) = i/t and r/(t) = r]/t. □ 
From Assumption 14.11 and the discussion from previous section, it follows that the 
optimizer of each maximization is unique and denoted by 

Xi{t, A) := argmax{-/i(xi) - tcpxA^i) - i^^ ^i^i) ■ ^ int(^i)) ^i^i = (14) 
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and x(t,X) := [xi(t, X)'^ ■ ■ ■ XNit, X)^]'^ . It is clear that the augmented dual function 
(P'^{t, X) = —d{t, X) and let X{t) := argmaxAgR™ (P'^{t, A), or equivalently 



X(t) = arg min d{t, X) 



From Assumption 14.11 and the proof of Lemma 13.31 it follows that the Hessian 
V'^dit, X) is positive definite for all t > and A G M™. Hence, the dual function 
d{t, ■) is strictly convex and thus X(t) is unique. Therefore, we can consistently define 
the set {{x{t, X{t)), X{t)) : t > 0}, called the central path. Let us introduce the iV^-self- 



concordant barrier function (px{x) '■= YliLi 4>Xi{xi) for the set X, where = XliLi 



Lemma 4.2. The central path {{x{t, X{t)), X{t)) : t > 0} converges to the optimal 
solution {x*, A*) as t ^ and {x{t, A(t)) : t > 0} is feasible for the problem 

Proof. Let x(t) := argmin2,.{/(x) + t(f)x{x) : Bx = b,Xi E int(Xj), AiXi = Cj Wi}, then it 
is known that x{t) — x* as t — 0. It is easy to see that the Hessian of f+t(j)x is positive 
definite and thus / + t(j)x is strictly convex and x{t) is unique. From Assumption 14.11 
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it also follows that strong duality holds for this barrier function problem and therefore 



min{/(a;) + t(px{x) : Bx = b,Xi E int(Xj), AiXi = ai Vz} = 

X 

maxmin{/(x) + t(()x{x) + (A, Bx — b) : Xi E int(Xj), AiXi = ai Vz} = 

A X 

min{/(x) + t(j)x{x) + (A(t), Bx — b) : Xi E int(Xj), AiXi = ai Vi}. 

X 



In conclusion, x{t) = x{t,\(t)) and thus x(t,X{t)) x* as t 0. As a consequence 



it follows that x{t,X{t)) is feasible for the original problem, i.e. Bx{t,X{t)) = b, 



AiXi{t, X(t)) = tti and X{t)) G int(Xj). It is also clear that X{t) — t- A* as t — t- 0. □ 



The next theorem describes the behavior of the central path: 



Theorem 4.1. For x{t) = x{t,X(t)) the following bound holds for the central path: 



given any < t < t then, 



f(x{t))-f{x{T))<N^{t-T). 
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Proof. For any s > 0, x{s) = [xi{s)'^ ■ ■ ■ xn{s)'^]'^ satisfies the following optimality 
conditions (see (|TT|) and Lemma [4 ■2p : there exists z/(s) G M^»=i™^ such that 



V/(x(s)) + sV4>x{x{s)) + B^\{t) + D\u{s) = 0, Bx{s) = b and AiXi{s) = ai. 



It follows immediately that (V f{x{s)),x'{s)) = —s{'V(f>x{x{s)),x'{s)). Since < 



T < t, then there exists t < s <t such that 



- /(x(r)) = (t-r)(V/(x(s)),x'(s)) 



s(t-r)(V0x(x(s)),x'(s)). 



From (j6]) we have that 



(V0x(x(s)),x'(s)) < (iV^ x'(s)^V'0x(x(s))x'(s)) 



1/2 
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Using a similar reasoning as in Lemma 13.41 we have: 



x'{s) = -[H~\s) - H-\s)Dl{DAH~\s)Dly^DAH-\s)]W(Px{x{s)), 

where we denote with H{s) = V^/(x(s)) + s'V'^(j)x{x{s)). Using ([H]), the expression 
for x'{s) and since -< V'^(j)x{x{s)) ^ l/sH{s) and H~^{s) ^ l/s(V^</)x(a:(s)))~"^ we 
obtain: 

x'{sfV^<f)x{x{s))x'{s) < {l/s)V(j)x{x{s)fH~\s)V<j)x{x{s)) 

< (l/s')V0x(x(s))^(vVx(x(s)))"V0x(x(s)) < N^/s^. 

It follows immediately that f{x{t)) - /(x(r)) < iV<^(t - r). □ 
A simple consequence of Theorem 14.11 is that the following bounds on the approxi- 
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mation of the optimal value function /* hold: 



< f{x{t)) -r< tN^. 

Indeed, from Lemma [4.21 we know that {x{t,X(t)) : t > 0} is feasible for the original 
problem ([I])-([2]). Since x{t) = x{t, X{t)), it follows that f{x(t)) > f*. It remains to show 
the upper bound. However, taking the limit as r — >■ in Theorem 14 . 1 1 and using Lemma 
14.21 we obtain also the upper bound. This upper bound gives us a stopping criterion in 
the algorithm that we derive below: if e is the required accuracy for the approximation 
of /*, then for any tf < e/N^ we have that x(tj) is an e- approximation of the optimum, 
i.e. x{tf) is feasible for problem ([I])-® and f{x{tf)) — f{x*) < e. Although X{t) is the 
minimizer of the dual function d(t, ■) over M™, so various unconstrained minimization 
techniques (e.g. Newton, quasi-Newton and conjugate gradient methods) can be used 
to approximate X{t), our goal is to trace the central path {{x{t, X(t)), X(t)) : t > 0} 
utilizing Newton method for the self-concordant family {d{t, ■)}t>o- 
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It is easy to see that the gradient of the self- concordant function d{t, ■) is given by 



TV N 



Vd{t, X) = b + Y^ Vdi{t, X)=b-Y^ BiXi{t, X) = b- Bx{t, A). 



i=l i=l 



For every (t, A) let us define the positive definite matrix 



H,{t, A) := V^fi{x,{t, A)) + tV''<PxM{t, A)). 



The Hessian of function di(t, ■) is positive definite and from (fT2l) it has the form 



W^d^it, A) = B,[Hiit, A)-i - H,it, Xy'Aj{AH,it, Xy'Aj) ' A,Hi{t, X)-']Bf. 



In conclusion, the Hessian of dual function d(t, ■) is also positive definite and given 

by: 

N 

V^dit,X) = ^V'rfi(t,A). 



i=l 
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Denote the Newton direction associated to self-concordant function d{t, ■) at A with 



AA(t, A) := -{V^d{t, A)) Vrf(t, A). 
For every t > 0, we define the Newton decrement of the function d{t, ■) at A as: 

S{t, A) := a{t)/2^Vd{t, X)T{V^d{t, A))^Vrf(t, A). 

Note that 6(t, A) = if and only if A = A(t) (recall that \{t) = argmiuAeR'" d{t, A)). 

Algorithm 4.1. (Initialization of Path- Following Algorithm) 
Step 0. input to > 0, Aq G M"", ey > and r = 

Step 1. compute x[ = Xj(to, A^) Vi, Sr = S(tQ, A^); if Sr < ey, rj = r and go to Step 3 
Step 2. determine a step size a and compute Newton iterate: Ar-+i = A^ -l- aAX{to, A^); 

replace r hj r + 1 and go to Step 1 
Step 3. output (t°,A°) = ito,Xrf). 
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Note that Algorithm 14. II approximates the optimal Lagrange multiplier A(to) of the 
dual function d{to, ■), i.e. the sequence (to, K) moves into the neighborhood V{t, ey) = 
{{t, A) : 6{t, A) < ev} of the trajectory {(t, A(t)) : t > 0}. 

Algorithm 4.2. (Path-Following Algorithm) 

Step 0. input: (t°, A°) satisfying 5(t°, A°) < ey , = 0, < r < 1 and e > 
Step 1. if t'^N^ < e, then kf = k and go to Step 5 

Step 2. ( outer iteration ) let = rt'^ and go to inner iteration (Step 3) 
Step 3. ( inner iteration ) initialize A = A^, t = t^"*"^ and 5 = 5{t''~^^, A'^) 
while 6 > ev do 

Step 3.1 compute Xi = Xi(t, A) Wi, determine a step size a and compute 
A+ = A + aAA(t,A) 

Step 3.2 compute 6^ = 6{t, A'^) and update A = A''" and 6 = 6^^ 
Step 4. A^"*"^ = A and x^~^^ = Xi] replace by + 1 and go to Step 1 
Step 5. output: (x;^^, ■ ■ ■ , , 
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In Algorithm 14. 21 we trace numerically the trajectory {{t, X{t)) : t > 0} from a given 



initial point A") close to this trajectory. The sequence {(x^, ■ ■ ■ ,x%, A'^)}a;>o lies in 



a neighborhood of the central path and each limit point of this sequence is primal-dual 



optimal. Indeed, since t^~^^ = rt^ with r < 1, it follows that lim,fc_^oo^'^ = and using 



Theorem 14.11 the convergence of the sequence = [{x^)'^ ■ ■ ■ (x^)'^]'^ to x* is obvious. 



The step size a in the previous algorithms is defined by some line search rule. 



There are many strategies for choosing r. Usually, r can be chosen independent of 



the problem (long step methods), e.g. r = 0.5, or depends on the problem (short 



step methods). The choice for r is crucial for the performance of the algorithm. An 



example is that in practice long step interior-point algorithms are more efficient than 



short step interior-point algorithms. However, short step type algorithms have better 



worst-case complexity iteration bounds than long step algorithms. In the sequel we 



derive a theoretical strategy to update the barrier parameter r which follows from the 



theory described in [16] and consequently we obtain complexity bounds for short step 



updates. Complexity iteration bounds for long step updates can also be derived using 
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the same theory (see Section 3.2.6 in [IE]). The next lemma estimates the reduction of 



the dual function at each iteration. 



Lemma 4.3. For any t > and A G M™', let AX = AX{t,X) be the Newton direction 



as defined above. Let also S = S{t, X) be the Newton decrement and = 2 — y/S. 



(i) If 6 > 5*, then defining the step length a = 1/(1 + 5) and the Newton iterate 



A"*^ = A + crAA we have the following decrease in the objective function d(t, ■) 



d{t, A+) - d{t, A) < -{4t/a^){6 - log(l + 6)). 



(a) If S < 5^, then defining the Newton iterate A"*" = A + A A we have 



S{t, A+) < 5V(1 - < V2, d{t, A) - d{t, X{t)) < {m/a^)6. 
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(in) If S < (5*/2, then defining = i^p[t, where c = 1/4 + 2^/6^, + rj, we have 

6{t+,X) < 6,. 

Proof, (i) and (ii) follow from Theorem 2.2.3 in [TB] and Lemma [4.11 from above. 

(iii) is based on the result of Theorem 3.1.1 in [16]. In order to apply this theorem, 
we first write the metric defined by (3.1.4) in [16] for our problem: given < < t 
and using Lemma 14.11 we obtain 

P5Mt^t+) = (l/4 + 2e/5. + r/)log(t/t+). 

Since 6 < 5^/2 < and since for = ^^^t, where c is defined above, one can 
verify that ps,/2it,t~^) = clog(l + l/2c) < 1/2 < 1 — 5/5*, i.e. satisfies the condition 
(3.1.5) of Theorem 3.1.1 in [T6], it follows that 6{t+, A) < 6,. □ 

Define the following step size: a{6) = 1/ (1 + 5) if 5 > 5* and cr(5) = 1 if 5 < 5*. With 
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Algorithm 14. 1 1 for a given t° and ey = 5*/2, we can find A°) satisfying 5(jp , A°) < 5*/2 
using tlie step size a{5) (see previous lemma). Based on the analysis given in Lemma 
143] it follows that taking in Algorithm |4i2] ey = ^*/2 and r = 2c/(2c + 1), then the 
inner iteration stage (step 3) reduces to only one iteration: 
Step 3. compute X''+^ = A^' + AA(t^+\ A^'). 

However, the number of outer iterations is larger than in the case of long step 
algorithms. 

4.2 Practical Implementation 

In this section we discuss the practical implementation of our algorithm and we give 
some estimates of the complexity for it. Among the assumptions considered until now 
in the paper the most stringent one seems to be the one requiring to solve exactly the 
maximization problems (1141) . i.e. the exact computation of the maximizers Xi{t, A). Note 
that the gradient and the Hessian of d{t, ■) at A depends on Xi{t, A)'s. When Xi{t, A)'s 
are computed approximately, the expressions for the gradient and Hessian derived in 
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the previous section for d{t, ■) at A are not the true gradient and Hessian of d{t, ■) at 



this point. In simulations we considered the following criterion: find Xi{t, A) G int(Xj) 



and z/i(t, A) G M™* such that AiXi(t, A) = and the following condition holds 



WVfiixiit, A)) + tV0x.(x^(t, A)) + BfX + Ajhit, A)|| < tc 



for some €3; > 0. Note however that even when such approximations are considered, the 



vector AA still defines a search direction in the A-space. Moreover, the cost of computing 



an extremely accurate maximizer of (fT^ as compared to the cost of computing a good 



maximizer of f[T^ is only marginally more, i.e. a few Newton steps at most (due to 



quadratic convergence of the Newton method close to the solution). Therefore, it is not 



unreasonable to assume even exact computations in the proposed algorithms. 
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4.2.1 Parallel Computation 



In the rest of this section we discuss the complexity of our method and parallel im- 



plementations for solving efficiently the Newton direction AA. At each iteration of the 



algorithms we need to solve basically a linear system of the following form: 



where d = Bi[Hr^ - Hr^Aj{AiHr^Aj) AiHr^]Bf, the positive definite matrix Hi 



denotes the Hessian of fi + t(j)Xi and some appropriate vector g. In order to obtain 



the matrices Hi we can solve in parallel N small convex optimization problems of 



the form (fT4|) by Newton method, each one of dimension Ui and with self-concordant 



objective function. The cost to solve each subproblem f|T^ by Newton method is 



0{n^{nx + loglog 1/tea;)), where nx denotes the number of Newton iterations before 



the iterates Xi reaches the quadratic convergence region (it depends on the update 



A) and te^ is the required accuracy for the approximation of (HM . Note that using 



n 




(15) 



i=l 
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the Newton method for solving (fT4l) we automatically obtain also the expression for 
and AiH^^Af. Assuming that a Cholesky factorization for AiH^^Aj is used to 



solve the Newton system corresponding to the optimization subproblem (fT4|) . then this 



factorization can also be used to compute in parallel the matrix of the linear system 



f|T5l) . Finally, we can use a Cholesky factorization of this matrix and then forward 



and backward substitution to obtain the Newton direction AA. In conclusion, we can 



compute the Newton direction AA in 0{^^^^nf) arithmetic operations. 



Note however that in many applications the matrices H^, Ai and Bi are very sparse 



and have special structures. For example in network optimization (see Section [S^ below 



for more details) the Hi's are diagonal matrices, Bi's are the identity matrices and the 



matrices Aj's are the same for all i (see (fT7|l ). i.e. Ai = A. In this case the Cholesky 



factorization of AH^ A can be done very efficiently since the sparsity pattern of those 



matrices is the same in all iterations and coincides with the sparsity pattern of AA^ , 



so the analyse phase has to be done only once, before optimization. 



For large problem instances we can also solve the linear system (1151) approximately 
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using a preconditioned conjugate gradient algorithm. There are different techniques 



to construct a good preconditioner and they are spread across optimization hterature. 



Detailed simulations for the method proposed in this paper and comparison of different 



techniques to solve the Newton system fll5l) will be given elsewhere. 



Let us also note that the number of Newton iterations performed in Algorithm 14.11 



can be determined via Lemma H^ (i). Moreover, if in Algorithm 14 . 2 1 we choose ey = 



and r = 2c/ {2c + 1) we need only one Newton iteration at the inner stage. It follows 



that for this particular choice for ey and r the total number of Newton iterations of the 



algorithm is given by the number of outer iterations, i.e. the algorithm terminates in 



polynomial-time, within 1^777^ ^og{N^t'^ / e)^ iterations. This choice is made only for 



a worst-case complexity analysis. In a practical implementation one may choose larger 



values using heuristic considerations. 
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5 Applications with Separable Structure 



In this section we briefly discuss some of tlie applications to which our method can be 



applied: distributed model predictive control and network optimization. Note that for 



these applications our Assumption 14. II holds. 



5.1 Distributed Model Predictive Control 



A first apphcation that we will discuss here is the control of large-scale systems with 



interacting subsystem dynamics. A distributed model predictive control (MPC) frame- 



work is appealing in this context since this framework allows us to design local subsystem- 



base controllers that take care of the interactions between different subsystems and 



physical constraints. We assume that the overall system model can be decomposed into 



appropriate subsystem models: 



x\k + 1) = ^ AijX^{k) + BijU^{k) \fi = l---N, 

jeAfii) 
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where A/'(z) denotes the set of subsystems that interact with the ith subsystem, including 



itself. The control and state sequence must satisfy local constraints: x^{k) G Qi and 



u^{k) G Ui for all i and k > 0, where the sets Qi and Ui are usually convex compact sets 



with the origin in their interior (in general box constraints). Performance is expressed 



via a stage cost, which we assume to have the following form: ii{x\ u^) , where 



usually £i is a convex quadratic function, but not strictly convex in Let Np 



denote the prediction horizon. In MPC we must solve at each step k, given x\k) = x\ 



an optimal control problem of the following form |18] : 



rnm[S^'S2^ii{x\,u\) Xq = x\ x\^^= V] Aijxj + Bijuj, x] e Qi, u] e Uiyi,i]. (16) 
1=0 i=i 



A similar formulation of distributed MPC for coupled linear subsystems with decou- 



pled costs was given in [3], but without state constraints. In [3], the authors proposed 



to solve the optimization problem (flGjl in a decentralized fashion, using the Jacobi al- 



gorithm [B]. But, there is no theoretical guarantee of the Jacobi algorithm about how 
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good the approximation to the optimum is after a number of iterations and moreover 



one needs strictly convex functions fi to prove asymptotic convergence to the optimum. 



Let us introduce Xj = [x^ - ■ ■ x^j^ Uq - ■ ■ u^j^_i\,Xi = fif^"*"^ x f//^ and the self-concordant 
functions /j(x*) = XlSo ^ (recall that ii are assumed to be convex quadratic). 



The control problem f|T6l) can be recast then as a separable convex program ([I])-(l2]), 



where the matrices Aj's and i?j's are defined appropriately, depending on the struc- 



ture of the matrices Aij and Bij. In conclusion, Assumption 14.11 holds for this control 



problem so that our decomposition method can be applied. 



5.2 Network Optimization 



Network optimization furnishes another area in which our algorithm leads to a new 



method of solution. The optimization problem for routing in telecommunication data 



networks has the following form [HE]: 



n N N 

{^fj{yj) + ^{ci,Xi):Axi = ai, '^Xi = y}, (17) 

j=l i=l i=l 
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where we consider a multicommodity flow model with commodities and n hnks. 



The matrix A G M"*^" is the node-hnk incidence matrix representing the network 



topology with entries { — 1,0, 1}. One of the most common cost functions used in the 



communication network literature is the total delay function (HE]: fj{yj) = ^ 



Corollary 5.1. Each function fj E C^{[0, dj)) is convex and fj is 3-compatible with 



the self- concordant barrier 4>j{yj) = — ].og{yj{dj — yj)) on the interval (0, dj). 



Proof. Note that the inequality ([9]) holds for all yj G (0, dj) and h eM.. Indeed, 



VV,(%)I = 3V'f,{y,)^l/{d,-y,y< f,{y,)^ 1/ {d, - y^Y + l/y^ □ 



Therefore, we can solve this network optimization problem with our method. Note 



that the standard dual function d^ is not differentiable since it is the sum of a differen- 



tiable function (corresponding to the variable y) and a polyhedral function (correspond- 



ing to the variable x). In [7] a bundle-type algorithm is developed for maximizing the 



non-smooth function do, in [T] the dual subgradient method is applied for maximizing 
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do, while in [6l[8] alternating direction methods were proposed. 



5.3 Preliminary Numerical Results 



We illustrate the efficiency of our method in Table 1 on a random set of problems of the 



form (fT7|) and (fT6|) . i.e. with total delay (ffist half table) and quadratic (second half) 



objective function, respectively. For the quadratic test problems we generate randomly 



the Hessian such that it is positive semidefinite of the form QfQi, where Qi are full 



row rank matrices. Here, the sets Xi are assumed to have the form [0, Ui] C i.e. 



rii = rii and also = mi for all i. Note that for these type of problems the barrier 



parameters Ni = 2ni and a < 8 and thus c = Ci + C2^/nl, for appropriate q > derived 



from Lemma 14.11 In our simulations we take r = 0.85, although a better tuning of 



this parameter will lead to less number of iterations. Complexity bounds for long step 



updates can also be derived using similar arguments as those given in the present paper 



for the short step method (see also Section 3.2.6 in [16]). For all test problems the 



coupling constraints have the form Xi = b, so that the total number of constraints 
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is equal to Nmi + rii. 



In the table we display the CPU time (seconds) and the number of calls of the dual 



function (i.e. the total number of outer and inner iterations) for our dual interior- 



point algorithm (DIP) and an algorithm based on alternating direction method [5] 



(ADI) for different values of mi,ni,N and fixed accuracy e = 10 ^. For two problems 



the ADI algorithm did not produce the result after running one day. All codes are 



implemented in Matlab version 7.1 on a Linux operating system for both methods. 



The computational time can be considerably reduced, e.g. by treating sparsity using 



more efficient techniques as explained in Section 14.21 and programming the algorithm 



in C. There are primal-dual interior-point methods that treat sparsity very efficiently 



but most of them specialized to block-angular linear programs [2]. For different data 



but with the same dimension and structure we observed that the number of iterations 



does not vary much. 
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Table 1: Computational results for network problems (fT7I) (first half) and quadratic 
problems ( IT6|) (second half) using DIP and ADI algorithms. 

6 Conclusions 



A new decomposition method in convex programming is developed in this paper using 
dual decomposition and interior-point framework. Our method combines the fast local 
convergence rates of the Newton method with the efficiency of structural optimization 
for solving separable convex programs. Although our algorithm resembles augmented 
Lagrangian methods, it differs both in the computational steps and in the choice of the 
parameters. Contrary to most augmented Lagrangian methods that use gradient based 
directions to update the Lagrange multipliers, our method uses Newton directions and 
thus the convergence rate of the proposed method is faster. The reason for this lies 
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in the fact that by adding self-concordant barrier terms to the standard Lagrangian 



we proved that under appropriate conditions the corresponding family of augmented 



dual functions is also self-concordant. Another appealing theoretical advantage of our 



interior-point Lagrangian decomposition method is that it is fully automatic, i.e. the 



parameters of the scheme are chosen as in the path-following methods, which are crucial 



for justifying its global convergence and polynomial-time complexity. 
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